Phase Transitions of S = 1 Spinor Condensates in an Optical Lattice 
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We study the phase diagram of spin-one polar condensates in a two dimensional optical lattice 
with magnetic anisotropy. We show that the topological binding of vorticity to nematic disclinations 
allows for a rich variety of phase transitions. These include Kosterlitz-Thouless-like transitions with 
a superfluid stiffness jump that can be experimentally tuned to take a continuous set of values, and 
a new cascaded Kosterlitz-Thouless transition, characterized by two divergent length scales. For 
higher integer spin bosons S, the thermal phase transition out of the planar polar phase is strongly 
affected by the parity of S. 



I. INTRODUCTION 

Low temperature superfluidity in simple atomic sys- 
tems is a well-studied subject. It is well-known that 
scalar bosons, confined to two spatial dimensions (2d), 
cannot form a true condensate. Instead, the off-diagonal 
correlations display either power law decay (in the super- 
fluid phase) or exponential decay (in the normal phase). 
These two phases are usually separated by a finite tem- 
perature Kosterlitz-Thouless (KT) driven by superfluid 
vortex unbinding. Such a transition in atomic systems 
has been observed experimentally through interference 
measurements.— 

Atomic systems of bosons with degenerate hyperfinc 
levels lead to more exotic phenomena at low tempera- 
tures usually referred to as "spinor condensation" . Such 
systems have been at the focus of intense experimental 
and theoretical activity since their discovery^ The hy- 
perfine levels give rise to a new quantum number anal- 
ogous to the spin. The macroscopic phase coherence in 
spinor systems can be accompanied by magnetic order. 
Indeed the spin and charge degrees of freedom may be 
strongly intertwined - as seen, for example, in the topo- 
logical defects, which can simultaneously involve atomic 
supercurrents and magnetic textures^ The presence of 
these multiple types of defects leads to a richer variety of 
phases and phase transitions. 

In this paper, we study thermal transitions in 2d polar 
condensates, where uniaxial spin nematic order (char- 
acterized by a headless vector, or "director" ) coexists 
with superfluidity First, we point out a crucial differ- 
ence between polar condensates with even integer and 
odd integer spin (as in S = 1 23 Na). For even S, the 
superfluid vortex and the nematic disclination are inde- 
pendent, while for odd S they are bound to each other 
topologically. This strongly impacts the phase diagrams. 
When the nematic director is confined to rotate in a 
plane, the normal state can be reached via a single con- 
tinuous transition from the polar state for the case of 
odd S, but not for even S, where a split transition is ex- 
pected. We study the odd S case in detail in this paper, 
and further show that the single continuous transition it- 



self can take on two very different characters, one, which 
is essentially KT-like, but with a non-universal superfluid 
stiffness jump; and another that is of a new 'cascaded- 
KT' type described in more detail below. Interference 
experiments- which have been used to study KT transi- 
tions in scalar condensates, can also be used to probe the 
new transitions discussed here. 

Consider spin-one bosons in a 2d optical lattice, de- 
scribed by a Hubbard model with couplings Uq and U2, 

H = -t 22 a i^ a 3°- + Uo /, n i (m — 1) — M /] n i 

(ij),cr i i 

+U 2 Y / (Sf-2n l )-gY / (SD 2 - W 
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The depth of the optical lattice serves to tune (Uo/t). 
Here a\ a creates an atom at site i with spin S? = a S 
{ — 1,0, 1}, rii is the particle number at site i, and [i is 
the chemical potential. The quadratic Zeeman field g, 
described below, is absent in magnetically isotropic sys- 
tems. We concentrate on atoms with antiferromagnetic 
spin interactions, U2 > 0, e.g. as is the case in 23 Na. 

The zero temperature phase diagram ofonodel |T]) with 
U2 > and g = was studied in Ref. [ij]- At unit fill- 
ing (one atom per lattice site), the system undergoes 
a continuous transition at T = between a nematic 
Mott insulator and a polar superfluid. The transition 
is tuned by the depth of the optical lattice. For deep 
lattices (Uo/t ^> 1) the system is a nematic Mott insula- 
tor, characterized by atoms which predominantly occupy 
the S z = state, together with a vanishing compressibil- 
ity. More generally, the nematic can be any h • S = 
state, with the director n serving as an order parameter. 
On the other hand, for weak optical lattices (Uo/t <C 1) 
the system is a polar superfluid. The order parame- 
ter = (a a ) = ("0+1, ipo, i ) -i) T m the polar state is 
* = e l6 TZ(0, 1, 0) T - % is a generic SO (3) spin rotation. 
As in the Mott insulator, the polar state has nematic 
order, described by a director n for which h ■ S\^) = 0, 
but it also has superfluid order, captured by an expecta- 
tion value of the superfluid phase 0. 

In the following, we are interested in systems with pos- 
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FIG. 1: Schematic phase diagram as a function of opti- 
cal lattice depth and temperature. The topological de- 
fects that disorder the polar state are: a superfluid vortex 
(q c , q a ) — (1,0), where 9 winds by 2n (red double circle); and 
a disclination+half- vortex (^, ^), where both 6 and cf> wind 
by 7T. Along the cascaded KT transition PQ, both defects 
play a role. 

itive quadratic Zeeman field g > 0. Such a field has the 
effect of restricting the director h to lie in the xy-plane 
in both the Mott nematic and polar states^ The most 
general planar polar state is then 

* = e ^(_ e ^ 7 O,e- i0 ) T , (2) 

where (j> is the angle of fi relative to the x-axis. The 
AC Zeeman effect - shining linearly polarized light at a 
frequency slightly detuned from the hyperfine level split- 
ting - can induce the required negative quadratic Zeeman 
field^ that leads to a planar polar state. The opposite 
g < limit, (induced by a magnetic field in S—l 23 Na) is 
essentially identical to a non-magnetic system since the 
nematic director is frozen along the field. 

II. TOPOLOGICAL DEFECTS 

Topological defects play an essential role in 2d finite 
temperature continuous phase transitions. In the present 
context, the planar polar SF is the "most ordered" phase, 
as it has both nematic and superfluid quasi long range 
order (QLRO). We can then understand the phase di- 
agram in terms of proliferation of defects of the planar 
polar condensate, which can destroy the order partially 
(leading to a nematic insulator) or completely. 

It is impossible to distinguish between the states ±n 
(n is a headless vector). However, an adiabatic rotation 
taking h to —h induces a change of sign in the polar order 
parameter This sign can be absorbed by simultane- 
ously shifting the superfluid phase 9 by ir. Therefore, to 
insure single-valuedness of the order parameter, in a po- 
lar state, a nematic disclination must be accompanied by 
a half-vortex in the superfluid phase. 

More generally, the topological point defects of a pla- 
nar polar condensate in 2d are labeled by two half- 
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FIG. 2: Topological defects of the planar polar state, (a) In 
a superfluid vortex (q c ,q 3 ) = (1,0), 6 C winds by 2-7T around 
the vortex (superflow shown as double circle with arrows). 

(b) In a nematic vortex (0, 1), the director h winds by 2tt. 

(c) In a disclination+half- vortex (A 1), h and 8 C wind by n 
simult aneously. 



integer charges, (q Cl q s ), describing the winding of 9 and 
</>, respectively, in units of 27r. By single-valuedness 
of \I> in Eq. the sum q c + q s is constrained to be 
an integer. The lowest energy defects are the super- 
fluid vortex (±1,0), the nematic vortex (0,±1), and 
the disclination+half- vortex (±|, ±i). These defects are 
shown in Fig. [5] 

Topological defects proliferate when the temperature 
is large relative to some appropriate stiffness parameter. 
In the case at hand, due to the presence of both spin and 
charge degrees of freedom, there are two relevant stiff- 
ness coefficients, K s and K c . These correspond to the 
energy cost of an elastic deformation in the nematic di- 
rection 4> and in the superfluid phase 9, respectively. In a 
dilute gas and in the absence of an optical lattice, the ki- 
netic term J^|W CT | 2 =|^|*| 2 {(V60 2 + (V0) 2 } predicts 
K c = K s . 6 On the other hand, enhanced quantum fluctu- 
ations in an optical lattice can change thisi An optical 
lattice suppresses both K c and K Sl but its main effect 
is to impede atomic motion, leading to K c /K s < 1. For 
strong lattice potentials at integer filling, i.e. in the Mott 
nematic phase, the charge stiffness is suppressed to the 
point where the system is an insulator, while maintain- 
ing nematic QLROi^ The nematic Mott insulator-polar 
SF transition is second order. Thus, proximity to this 
transition allows tuning the ratio K c /K s < 1 over a wide 
range. 

The topological defects interact logarithmically at long 
distances, leading to a Coulomb gas action, 

S = J2 (KoQtQj + K s qfq°) log %L + J2 log Vi (3) 

ij i 

Here, a is the defect core size, and the reduced stiffnesses 
K c = and K s = ^fr have been normalized by a factor 
of 2/7T for later convenience. The defect fugacity j/i takes 
the values y c and y m , respectively, for the defects (±1,0) 
and (±i, ±i) in Fig. H 



3 



III. DEFECT UNBINDING 

The phase diagram as a function of the reduced stiff- 
nesses K s and K c consists of four phases: (i) Polar 
state.- At large K s and K c , all defects remain bound 
and the polar order parameter has algebraic order, (ii) 
Disordered. - In the opposite limit, when both K s and K c 
are small, all of the topological defects proliferate, and 
the system has short-range correlations in both charge 
and spin, (in) Nematic- Starting in the polar state and 
keeping K s large, when K c is reduced sufficiently, super- 
fluid vortices proliferate, with all other defects remain- 
ing bound. This leads to algebraic order in the nematic 
order parameter e 2l< ^, but no superfluidity, (iv) Paired 
super fluid. - Conversely, when K c is large and K s is small, 
the nematic vortex unbinds, leading to algebraic order in 
e 2%e , but no spin order. Note that, for polar condensates 
in an optical lattice, for which K s > K c , the paired su- 
perfluid is not present, as shown in Fig.[TJ 

The phase boundaries of the polar state can be ob- 
tained from the requirement that all defects in Fig. [2] 
be bound. This corresponds to K*r > 1 (bound super- 
fluid vortices), Kg > 1 (bound nematic vortices), and 
K^ + Kg > 4 (bound disclination+half- vortices). Here 
we have introduced the notation Kf? s to denote the long- 
distance spin and charge stiffness, which is renormalized 
by the presence of a finite density of bound defect pairs. 
The conditions for vortex unbinding are shown in Fig. [3J 

To obtain the correct phase diagram, however, it is 
important to recognize that the conditions for defect un- 
binding cannot always be treated independently of one an- 
other. In particular, there are situations where the un- 
binding of one type of defect can precipitate the unbind- 
ing of a second type of defect. As an example, consider 
the triangular region enclosed by PQS in Fig. [31 At 
first glance, it looks like the superfluid vortices are pro- 
liferated (since K c < 1) whereas the disclination+half- 
vortices remain bound (since K c + K s > 4) . This would 
identify this region as belonging to the nematic phase. 
However, points within this region cannot describe a sta- 
ble nematic phase. The cheapest defects inside a nematic 
are the single disclinations - in which the phase <f> winds 
by 7r. Disclinations are bound whenever the reduced spin 
stiffness K s is larger than 4. However, K s < 4 in the re- 
gion PQS. Therefore, this region is part of the disordered 
phase. 

The instability of the nematic within region PQS can 
also be understood by thinking of the single disclination 
as the remnant of the disclination+half- vortex of the po- 
lar phase, once the superfluid phase has been disordered 
due to the proliferation of superfluid vortices. In this sit- 
uation, the proliferation of superfluid vortices renormal- 
izes the charge stiffness to zero (K^ = 0). Therefore, al- 
though the bare values of K c and K s satisfy the condition 
for disclination+half-vortex to be bound, K c + K s > 4, 
the renormalized stiffnesses do not, Kf -\-Kf < 4. These 
arguments show that, along the line PQ (and also along 
P'Q'), the system undergoes a direct transition from the 




FIG. 3: Phase diagram as a function of the reduced stiff- 
nesses Kc and K s . The solid lines represent the conditions 
for the different topological defects to be bound, when the 
interactions between different types of defects are ignored. 
Two regions (enclosed by PQS and P'Q'S') are labeled by 
question marks. In the absence of vortex interactions, these 
regions would be in the nematic and paired superfluid phases 
(pSF), respectively. However, as argued in the text, the or- 
dered phases in these two regions are unstable, and the two 
regions are part of the disordered phase, which extends to the 
dashed lines QS and Q'S'. For simplicity, the phase diagram 
is drawn for infinitesimal bare defect fugacities. Point O is 
not a special point in the phase diagram, but it is indicated 
for comparison with Fig. [T] 



polar state to the disordered state. This transition in- 
volves two different types of defects, and therefore is in 
a different universality class from the transition along 
P'OP. We call the transition along PQ a "cascaded 
Kosterlitz-Thouless" transition (cKT). 

It is difficult to give a fully controlled renormalization 
group (RG) treatment of the cKT transition, since the 
physics involves the proliferation of one type of defect be- 
fore the other defect "realizes" that it is unbound. The 
conventional RG treatment of defect unbinding is only 
controlled when the defect fugacities are small, and it 
breaks down when the defects proliferate and the fugac- 
ities become large. Note that, despite the difficulty in 
giving a carefully controlled RG treatment, the coarse 
graining process in the RG cannot increase the value of 
the stiffnesses. Therefore the argument for the instability 
of nematic order within the region PQS is robust. In the 
following section, we will study the RG equations and de- 
rive the critical properties of the cKT transition PQ, and 
we will contrast it to the disclination+half-vortex prolif- 
eration along P'OP. This will be followed in Sec.|V]by a 
numerical study of the transitions along PQ and P'OP 
using Monte Carlo simulations. This will allow us to 
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study the finite size scaling properties at the two transi- 
tions, and compare them with the predictions from the 
RG analysis. 



IV. RENORMALIZATION GROUP ANALYSIS 

A. Renormalization group equations 

In this section we focus on the case K s > K c , cor- 
responding to bosons in an optical lattice. Then, the 
nematic vortices (0, 1) are always the last defects to pro- 
liferate, and they can be neglected in the analysis below. 
A real space renormalization group (RG) analysis of the 
Coulomb gas ([3]) is carried out in Ref. [8(. To quadratic 
order in the fugacities, § 

y c =2(1 - K c )y c + 2ny 2 m (4a) 

y m =(4 - K s - K c )y m /2 + 2ny m y c (4b) 

k c =-8nK 2 c (2y 2 c +y 2 m ) (4c) 

K s =- ^K 2 s y 2 m (4d) 

Here g = % and e e is the length rescaling factor. These 
RG equations are valid provided the fugacities remain 
small. As is clear from Eq. ([3]), the RG equations must be 
symmetric under the exchange of spin and charge degrees 
of freedom. 

In order to understand the RG equations ([4]), consider 
first the flow of the superfluid vortex fugacity y c , Eq. (|4"aj) . 
The first term on the RHS describes the competition be- 
tween energy cost and entropy gain for the creation of a 
superfluid vortex; the second term, proportional to y^, 
arises because two dislocation+half-vortices can combine 
into a superfluid vortex. On the other hand, Eq. (|4"c)) for 
the flow of K c describes the screening of the Coulomb in- 
teraction between defects due to a finite density of bound 
defect-anti-defect pairs in the medium. The sign of the 
flow of K c is negative semi-definite, and is only zero when 
the fugacities y c and y m are zero, or when K c itself is zero. 
Thus, the fixed point in the RG requires either the fugaci- 
ties or K c to flow to zero. Similar considerations apply to 
K s . The fixed points of the RG are characterized by the 
values of the renormalized stiffness Kf = K 7 (£ = oo), 
(7 = c,s). 

The continuous transitions between the phases in 
Fig. [1] arise from defect unbinding, and may be classified 
according to the type of defect that triggers the transi- 
tion. If a single type of defect is important, one observes 
a conventional KT scenario. However, we also find a class 
of transitions where unbinding of one set of defects trig- 
gers the instability in another set, leading to a cascaded 
KT transition with two diverging length scales. Here, we 
will concentrate on the two direct transitions between the 
polar and disordered states, along the lines OP and PQ. 
The multicritical point P belongs to a different univer- 
sality class, as discussed in Ref. |8|. The other transitions 
QR and QS belong in the conventional KT scenario. 



B. Disclination+half-vortex unbinding (OP) 

Here, the polar state is disordered by the prolifera- 
tion of disclination+half-vortices. These defects destroy 
both charge and spin order. At the transition, the renor- 
malized stiffness satisfy Kf + Kf = 4 and Kf > 1. 
The superfluid stiffness jump in this transition can be 
tuned continuously with optical lattice depth, ranging 
from Kf = 2 at point O, to Kf — > 1 as we approach 
point P. On the other hand, the sum of superfluid and 
spin stiffness is universal. The correlation lengths £ 7 di- 
verge as the transition is approached from the disordered 
side as 



a exp 



d 7 (T 



Tkt) 



(5) 



As in the usual KT transition, d r = d. is non-universal. 



C. Cascaded KT criticality (PQ) 

Along PQ, the superfluid vortices are on the verge of 
proliferating, since Kf = 1. On the other hand, the 
sum Kf + Kf is above the threshold value of 4, indicat- 
ing that the disclination+half-vortices are bound at the 
transition. However, as soon as the transition is crossed, 
K c flows to zero, reducing Kf + Kf below 4. Now, 
the disclination+half-vortices unbind, leading to a com- 
pletely disordered phase. We call this a "cascaded" KT 
(cKT) transition, since unbinding of one type of defect 
triggers the unbinding of the other. 

Both spin-nematic and charge orders have diverging 
correlation lengths as PQ is approached from the dis- 
ordered phase. However, there is a separation of scales 
S,s 3> £0 due to the fact that the superfluid vortices un- 
bind at a shorter length scale than the disclination+half- 
vortices. We will show below that the two are related by 
a power law, 



afe/a) 1 



(6) 



where B = 1/(4- Kf) > 1, and £ c follows Eq. ©. 

The main challenge in studying the cascaded KT tran- 
sition is that the naive RG equations ^ break down 
once the superfluid vortices unbind. In order to circum- 
vent this, we use the separation of scales to perform the 
RG in two steps. First, we solve Eq. ((4]) up to the scale 
£ c where the superfluid vortex fugacity begins to diverge. 
At this point, the superfluid correlations are explicitly 
short-ranged, and we can integrate out the charge de- 
grees of freedom to obtain a local spin-only model. We 
then study the RG flow of the ensuing spin model. 

In the first step, the fugacity of the disclination+half- 
vortex is renormalized down according to Eq. (|4b[) . y m ~ 
ym{ a /£c)^ K "~ 3 " 2 - The RG flow of the coarse grained 
couplings y m , and K c at longer scales is then governed by 
Vm = \ (4 - K^j y m and K s = -8-Kk 2 s y 2 m . Integrating 
this RG flow until y m is of order unity yields Eq. ((5]), 
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FIG. 4: An illustration of a current loop configuration dis- 
cussed in the text. The bonds with arrows are the directed 
dimers and solid undirected bonds are the undirected dou- 
ble dimers. The configuration shown is planar, while in the 
model studied here the configurations will be similar but on 
a cubical lattice. 



with B = (4 - Kfy 1 . Note that B is bounded below 
by one, near point P, and can be arbitrarily large near 
Q, where Kf = 4. The topology of the phase diagram 
in Fig. [1] is crucially different from Ref. 0] , where the 
cKT transition (PQ) is misinterpreted as two separate 
transitions. 

The two diverging length scales at the cKT transition 
arise due to the fact that the disclination+half-vortex is 
a dangerously irrelevant operator at the transition. This 
is one of the few examples we know of a dangerously 
irrelevant disorder operator j&iii 



V. MONTE CARLO SIMULATIONS 



A. Loop model 



Roughly, the parameter t is a temperature- like parameter 
(raising t eliminates in-plane loops and double dimers), 
while Wo tunes the ratio K c /K s (increasing Wo in- 
creases the double dimer density). As an illustration we 
show a planar lattice configuration in figure[4j Configura- 
tions that contribute to the partition function are similar 
configurations but on a cubic lattice of size L x L x 4. 

The loop model ((7|) is a simple variant of strongly 
coupled two color lattice QCD with staggered fermions 
(SCLQCD2), a model with an 50(3) x U(l) symmetry 
and hence of interest also in the field of spinor conden- 
sates. In the model we consider here, the 5*0(3) symme- 
try is broken to an 50(2) subgroup. While SCLQCD2 
has been studied in both cubic (3d)ii and hyper-cubic 
(4d)^2 lattices, the above model remains unstudied so far. 
These models can be studied efficiently using directed- 
loop Monte Carlo algorithms as discussed in Refsfi~2l 
andO 

As a consequence of the constraints, the model con- 
sidered here has two conserved currents: (1) Jf — 
Vi{\ e i,a\ +^i,q ~ 1/3} where r/i = +1 and —1 on alternat- 
ing sites, and (2) J? a = e^ a . These currents correspond 
to the spin and charge conservation and can be used to 
compute K s and K c : 



K- 



(8) 



Here, K 1 is normalized such that at a usual KT transi- 
tion one would expect K s , K c = 1 . Note that at every 
space-like link, Jf a + J? a is an even integer—. This im- 
plies that the topological defects are half-integers with 
the constraint that q c + q s is an integer, which is cru- 
cial to the physics hereSi. In Appendix lAl we provide a 
diquark representation of the loop model, in which the 
symmetries of the model are seen explicitly. Due to these 
symmetries, the model studied here is expected to exhibit 
the same universal physics as the Coulomb gas system de- 
scribed by Eq. §5§ close to second order phase transitions. 
In fact we will provide below clear numerical evidence 
for the predicted transitions along the line OP (discli- 
nation+half vortex unbinding) and PQ (cKT transition) 
using this model. 



In order to test the new transitions predicted in the 
previous sections, we introduce a relatively simple lattice 
current-loop model defined on a periodic cubic lattice. 
The configurations [C] that contribute to its partition 
function contain two types of bond variables: directed 
dimers e.; jQ , € {1,0,-1} and undirected double dimers 
di.a € {0, 2}, on bonds between site i and its neighbors i+ 
a (a — x, y, r). At each site i we enforce two constraints: 
(1) Eo{ e i," + ei-&,a\ = (directed dimer conservation) 
and (2) J2 a {\ e i,a\ + d t . a + \ei^ &<ol \ + d l - & . a } = 2 (close 
packing constraint). The partition function is given by 

Z = II {Wo} {< * +di ' v+d ^ )/2 {t}l ei -l+*- . (7) 
c i 



B. Disclination+half-vortex proliferation (P'OP) 

Let us first focus on the transition that occurs along 
the line P'OP. When one is exactly on this line, the finite 
size scaling formula for K a , can be computed using RG 
and is given by 

where + Kf — 4 and the C a (a = c, s) are non- 
universal constants. We would like to show evidence for 
this using our model. We have discovered that in our 




FIG. 5: Finite size scaling near the polar melting transition (P'OP) at Wd = and t = 0.9, 1.0, 1.05, 1.1 and 1.2. The solid 
line is the fit to Eq. ©. The values of the fit are given in table U 



model we can approach this line by varying t for fixed 
Wjj in the interval < Wn < 3. We have done extensive 
calculations up to lattice sizes of L = 512 by focusing on 
Wd = and 1. In Fig. [5] we show the finite size scaling 
of K s and K c at Wd — for various values of t close 
to the transition. Tabic Q] shows that we obtain a good 
fit to Eq. @ for t < 1.05. This is understandable since 
everywhere inside the polar phase K a is expected to be 
a constant for large L and this just means C a is large. 
Since t = 1.05 is the last value of t where the fit works 
well, at that value we should be very close to the critical 
line. It is important to note that indeed K c + K s 4 
there (see last column of table Q1 . Individually each of 
the ^'s are not universal. Since Kf > Kf we must be 



t 


K? 


C c 


x 2 


Kf 


c s 


x 2 


K* + K* 


W D = 


0.80 


2.96(2) 


7(1) 


0.3 


1.95(2) 


oo 


0.7 


4.9(4) 


0.90 


2.74(1) 


6.5(9) 


0.3 


1.89(1) 


10(2) 


0.9 


4.63(2) 


1.00 


2.468(5) 


3.1(2) 


1.6 


1.763(5) 


4.7(5) 


1.1 


4.231(10) 


1.05 


2.297(2) 


1.09(3) 


0.9 


1.677(2) 


2.05(7) 


0.7 


3.974(4) 


1.10 


2.145(2) 


0.34(2) 


24.7 


1.594(2) 


0.98(4) 


4.4 


3.739(4) 


Wd = 1 


0.80 


2.651(1) 


oo 


1.4 


2.652(1) 


oo 


1.3 


5.316(3) 


0.90 


2.41(2) 


15(5) 


0.6 


2.40(1) 


12(3) 


1.1 


4.81(3) 


1.00 


2.138(6) 


4.2(4) 


1.3 


2.144(9) 


4.6(7) 


2.3 


4.282(15) 


1.05 


1.967(2) 


1.29(3) 


2.6 


1.968(2) 


1.29(5) 


2.3 


3.935(4) 


1.10 


1.821(1) 


0.43(2) 


59 


1.818(2) 


0.40(2) 


24 


3.639(3) 



TABLE I: Fits to Eq. % at Wd = and Wd = 1. Note that 
t = 1.05 is almost on the P'OP line, with + Kc » 4 with 
goodness of fit x 2 /DOF wl-2. When W D = 1 and t = 1.05 
we are at the midpoint O on the P'OP line. 



somewhere on the P'O line of the phase diagram. We 
can change WD and t to go to a different point on the 
P'OP line. To show this have performed calculations at 
WD = 1. In table [J we also give some of our fits for the 
WD = 1 case. Again t = 1.05 seems to on the critical 
line and K c + K s w 4. But now we have K c = K s which 
means we are right at the midpoint O. 



C. Cascaded KT transition (PQ) 

Next we focus on the cascaded KT transition line PQ. 
On this line the cKT scenario predicts that K s and K c 
should follow different finite size scaling laws: 



K c (t cKT ,L) = l + i__L 

2 C c + log L 

K s (t cKT ,L) = K? + C S L 3 - K ?, 



(10) 
(11) 



where C 7 are again non- universal constants, and 3 < 
< 4. In order to test this scenario, the first step 
is to locate a point on the cKT line within our model. 
We have identified numerically that one such point is 
Wd = 3.05 and t — t C KT = 0.885. In order to cross the 
transition line at this point we vary both Wd and t such 
that W D = 3.05 + 10(t- 0.885), and compute K s and K c 
for lattice sizes up to L = 512. 

Figure [5] shows our data for both K c and K s . The 
first important qualitative observation we make from the 
figure is that both K s and K c appear to jump to zero 
for t > t C KT as expected. Further the rough value of 
the jump is rs 1 and Kf" w 3.6 which is again as 
expected. Thus, the point is roughly midway between P 
and Q, yet in the numerical data, both K s and K c seem 
to undergo a transition simultaneously. This is consistent 
with the results of Sec. lIIIi where we argued that for these 




FIG. 6: Cascaded KT (PQ): K c (top) and K„ (bottom) vs. t/t C KT for various system sizes L. Insets: K c and K s vs. L at 
fixed t = t C KT, fit to Eqs. (JTUJ) and flTTJ, respectively. The fits yield Kf = 3.567(3), C s = 0.58(2), and C c = -0.168(4), with 
goodness of fit x 2 /DOF ~ 1. 



parameters a split transition i.e. crossing QR first and 
then QS is ruled out since Kf < 4 is inconsistent with a 
stable nematic phase. 

Next we perform a more quantitative analysis at t = 
t C KT- We fit our data for K c and K s as a function of 
L to Eqs. ([TDJ and JTIJ. If iff = 1 is fixed we hnd 
C c = -0.168(4), with goodness of fit X 2 /DOF = 1, and 
Kf = 3.567(3), and C s = 0.58(2) with X 2 /DOF = 0.7. 
We emphasize that 8 < L < 512 were used in the fit. 
This large range gives us confidence in our analysis. The 
insets in the graphs in Fig. [6] show the finite size scaling 
of K c and K s , respectively along with the fit. Thus, we 
claim that the Monte Carlo simulations are consistent 
with the expected finite size scaling predicted by the RG 
analysis. 

Another important prediction of the cKT scenario is 
the presence of two diverging correlation lengths near the 
cKT transition PQ, which follows Eq. ©. Unfortunately 
we have not measured the correlation lengths directly. 
However we can see this divergence in an indirect manner. 
Since a finite size box limits the diverging correlation 
lengths, the jumps in the stiffness K a are no longer sharp 
but broadened depending on the value of L (see Fig. [5]) . 
A temperature t* > t^r can then be defined such that 
K a {t* a ) — xKf(tKT) where x < 1 is some fixed fraction. 
In other words, at t*, K a has been reduced to a fraction 
x of its value at txr- The fact that the correlation length 
diverges as t approaches tjcr is now seen by the fact that 
as L becomes large t* approaches txT- This behavior 
can be quantified using RG and we expect 



1 



a a log L + p a 



(12) 



• 


t c "(x=0.5) 


o 


t c '(x=0.1) 


A 


t s *(x=0.5) 


□ 


t 8 *(x=0.1) 




where all constants are non-universal, except for the ratio 
a c /a s = l/(A-Kf). Note that Eq. flT2j! is valid provided 



FIG. 7: (t* -iA'r)~ 1/2 vs. system size for charge correlations 
and spin correlations at x — 0.1 and x = 0.5. The solid lines 
are fits given in table Ull 



that i* — tKT is small, i.e. for large enough system size 
L. 

Fig. [7] shows our data for t* for the choices x = 0.1 
and x = 0.5. The first striking qualitative observation 
we make is that 1/y/t* — Tkt scales linearly with log(L) 
to a very good approximation as expected. Further the 
slopes of the two lines are clearly different. Note that 
we expect £ s ~ £^ where p — a c /a s > 1. Indeed we 
do find that p > 1. Quantitatively, while t* fits well to 
Eq. (fl"2")) . t* does not fit for larger values of L due to a 
small but clear curvature in the data. Table |TT] contains 
the fit results. Ignoring the large values of % 2 /DOF, we 
see that a c /a s ~ 3.5 at x = 0.1 and ~ 3.2 at x = 0.5. 
These must be compared to the theoretical prediction 
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X 


a„ 




X 2 /DOF 


ttc 




X 2 /DOF 


0.1 


0.321(2) 


0.577(8) 


8.3 


1.125(5) 


0.20(2) 


1.4 


0.5 


0.363(3) 


0.74(1) 


35 


1.15(1) 


1.25(4) 


0.3 



TABLE II: Fits to Eq. at x = 0.1 and x = 0.5. The fits 
for ts are not very good because of a detectable curvature in 
the data. 



from the RG treatment, 1/(4 - Kf) = 2.31(1). If we 
take the curvature in t* into account by using a fit of the 



form 



= a s logi + /3 s + 7s /(logL) the x 2 /DOF 



improves slightly but is still not good. The ratio a c /a s 
changes to about 3 and 2.3 at x = 0.1 and 0.5. If in- 
stead we just use the last two points and draw a straight 
line through the data, this ratio is 2.8 and 1.85 respec- 
tively. This large variation in the values of the ratio 
a c /a s shows that we do not yet have quantitative con- 
trol on it. The theoretical expectation seems to be within 
the large systematic errors. Since our data are obtained 
on rather small lattice sizes where the two large length 
scales £ s ~ £ 2 ' 3 cannot fit well, our inability to quantita- 
tively control the ratio is not surprising. 



VI. CONCLUSION 

In conclusion, we have shown that the topological bind- 
ing of spin and charge vorticity in S = 1 polar conden- 
sates can give new types of phase transitions. In particu- 
lar we have shown the existence of the disclination+half 
vortex unbinding transition, which is similar to the KT 
transition but where the superfluid stiffness jump is non- 
universal and a cascaded KT transition where in addition 
to a non-universal jump in the stiffness two correlation 
lengths diverge at the critical point, where one is a power 
of the other. A large scale numerical study supports the 
detailed picture we have developed. 

Our analysis can be applied to other ordered 
states. For example, spin-2 87 Rb in a magnetic 
field forms a "square nematic" stat o 15 > 16 with W ~ 



e ia (e 2l f>, 0,0,0, 



! ). This yields precisely the same 



physics as the spin-1 planar polar state. There are exam- 
ples outside of cold atomic systems which display the the 
same topological structure studied here. For instance, 
Berg el al. have argued that the thermal melting of a 
striped superconducting state is produced by a fully anal- 
ogous set of topological defects to those discussed here. 17 
Therefore, the phase diagram of a striped superconductor 
is the same as that of a planar polar condensate, although 
the particular phases involved in the two cases are differ- 
ent. Following this work, a recent numerical study has 
looked at the 3d version of this system*^ 

We thank E. Berg, G. Delfino, J.E. Moore, S. Muker- 
jee, L. Radzihovsky, and D. Stamper-Kurn for useful dis- 
cussions. This work was supported in part by the NSF 
grant DMR-0506953 and the Hellman Faculty Fund. 



APPENDIX A: DIQUARK REPRESENTATION 
OF THE LOOP MODEL 



In this Appendix, we recast the loop model Eq. |J7J| in 
terms of a path integral over Grassmann variables in or- 
der to clarify the symmetries of the model. The partition 
function is given by 



(Al) 



where the action 



Sc 



{(* i+a ^) 2 + (^*i +a ) 2 } 



Here 



9i 




and $>i 



it til 



(A2) 



such that ipia and i>ia {°~ = t>4) are independent Grass- 
mann variables residing on the sites of a L x L x 4 lattice. 
In particular, ipia and ip ia . are not complex conjugates of 
each other. The constants 1 and 

u T = t, so that f can be thought of as an Euclidean time 
direction and t as a temperature-like parameter. Per- 
forming the path integral over tp^ and tp ia yields the 
loop model ([7]) exactly, as can be checked easily by ex- 
plicit computation. 

Equation (jAlj) can also be expressed as a model of 
"diquarks" Di = ipiii/j^ and Di = ^tV'ii hopping and 
interacting on a lattice, 

Sg = )]u a (AA+a + A+aA) 

- (W D - l)J2 u lDiDiD i+a D i+a . (A3) 



The model has an SU(2)xSU(2) gauge symmetry. To see 
this, note that the diquark Di — ipiiipi] — ^aa'^ia^ia 1 
is invariant under a local transformation 



(A4) 



where Ui G SU(2). Similarly, an independent SU(2) 
transformation Ui can be carried out on the barred vari- 
ables, ^i — > \&j u\ 7 leaving Di invariant. 

After taking the gauge invariance into account, we are 
still left with an independent global U(l)xU(l) symmetry 
of the model (|A1[) , parametrized by the angles and <p, 



(A5) 



9 



where rji — +1 and —1 on alternating sites. Note that, 
unlike the gauge invariance in Eq. (|A4|1 . the transfor- 
mation (|A5|) does not leave the diquark invariant, since 
V'iJ.V'iT ~~ * ^^^'^ipiiipn; . However, a simultaneous shift 
of 9 and <f> by the angle n does leave the diquark invari- 
ant. Hence, the model (|Alj) is explicitly seen to have 
topological defects that are labeled by half-integers with 



the constraint that q c + q s is an integer. Therefore, the 
low energy defects of the model are the superfluid vortex 
(1,0), the nematic vortex (0,1), and the nematic+half- 
vortex (i, i). By equivalence of the diquark model to 
Eq. ([7]), we conclude that the loop model also has the 
same topological defects. 
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